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ABSTRACT 

A systematic procedure for constructing semidiscrete, second order 
accurate, variation diminishing, five point band width, approximations to 
scalar conservation laws, is presented. These schemes are constructed to also 
satisfy a single discrete entropy Inequality. Thus, in the convex flux case, 
we prove convergence to the unique physically correct solution. For 
hyperbolic systems of conservation laws, we formally use this construction to 
extend the first author's first order accurate scheme, and show (under some 
minor technical hypotheses) that limit solutions Satisfy an entropy 
inequality. Results concerning discrete shocks, a maximum principle, and 
maximal order of accuracy are obtained. Numerical applications are also 
presented . 
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0. Introduction 


Recently a number of new shock capturing finite difference approximations 
have been constructed and found to be very useful in shock calcu]a,tions, e.g. 
[3]j [^]j [16]. In addition to conservation form, these schemes are usually 
constructed to have as many as possible of the following properties: 

(1) Stable and sharp discrete shock solutions 

(2) Limit solutions which satisfy a geometric and/or analytic entropy 
condition. 

(3) A bound on the variation of the approximate solutions, at least in 
both the scalar, and linear systems, case. 

(h) Second order accuracy in regions of smoothness, (with certain 
isolated exceptional points, as described in Section II below). 

In this pajier we shall present a general procedure for constructing schemes 
with a five point band width satisfying all of the above. We shall then 
prove a convergence theorem for a wide class of approximations to any scalar 
convex conservation laws in one space dimension, using recent uniqueness 
results of DiPerna [5]. 

Our convergence proof Involves the following simple steps: 

(a) We first obtain a variation bound for a wide class of second order 
accurate approximations. This implies that, for any sequence of mesh widths 
approaching zero, a subsequence of approximate solutions converges to a weak 
solution of the Cauchy problem. 

(b) Next, we obtain a single discrete entropy inequality for a large 
sub-class of the approximations mentioned above. Then all limit solutions 
in this subclass satisfy the inequality. 
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(c) We finally invoke the uniqueness results of DiPerna to prove 
convergence to the solution to the Cauchy problem, as the mesh width goes to 
zero, for the subclass mentioned in (b). 

Earlier, Lfejda and Osher [20], modified the Lax-Wendroff scheme, keeping 
its essential properties, so that an entropy inequality was obtained. DiPerna 
(private communication) has proven IT convergence of this modification 
(again approximating convex conservation laws) using the theory of compensated 
compactness. No variation bound is possible for this second order accurate 
approximation. 

From a practical point of view, this lack of variation bound means that 
conventional schemes such as Lax-Wendroff, when approximating hyperbolic 
systems of conservation laws, even with an entropy modification, suffer from 
a lack of robustness in computing complex flows with shock waves and steep 
gradients. While such schemes have been widely used in a variety of problems, 
(see [2] for references), that list of solved problems does not include flows 
with strong shocks (say Ifech 5 upstream, normal to the shock), when the shocks 
are captured. 

A main drawback of most finite-difference schemes is that discontinuities 
are approximated by discrete transitions, that when narrow, usually overshoot 
or undershoot, or when monotone, usually spread the discontinuity over many 
grid points. 

Upwind schemes have been designed and used over the years, largely 
because of their success in treating this difficulty. Those based on solving 
the Riemann problem either exactly (Godunov's method [10]) or approximately 
e.g. (Osher's [2^], or Roe's {28] with an entropy fix [3^], [22]), have 
been extremely successful, especially when made second order accurate, [3], 
[4], or [ 16 ]. 
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We should particularly mention the early Investigations of van Leer [ 32 ], 
[33]. Hiere he introduced the concepts of flux limiters, and higher order 
Riemann solvers. Recently, Harten [I6], using an argument also used in 
[1], [31] and [29]) obtained simple sufficient conditions compatible with 
second order accuracy, which guarantee that a scalar one dimensional approxi- 
mation is TVD - total variation diminishing. He constructed a scheme having 
that property, and formally extended it to systems, using a field -by -field 
limiter, and Roe’ s decomposition. 

Using any of the three nonlinear decompositions (Godunov' s, Osher' s, 
or Roe's) one can obtain field-by-field limiters for systems, as described 
in Section VI below. Peter Sweby, [30], has investigated the properties 
of various limiters, and clarified their application considerably. The notation 
used in Sections III -VI below is due to him. 

We shall u£se the now introduced term "high resolution scheme" to mean 
a formal extension to systems via a fleld-by-field decomposition, of a scalar, 
second order accurate, variation diminishing scheme as in [16]. These schemes 
do not, in general, satisfy the entropy condition - e.g. expansion shocks 
exist as stable solutions of high resolution schemes based on Roe's (unmodified) 
scheme. In Section VI, we use Osher' s decomposition and certain limiters, 
to prove, that under technically reasonable hypotheses, limit solutions of 
certain high resolution schemes do satisfy the entropy condition for hyper- 
bolic systems of conservation laws. 

Fundamental to our work is a formula measuring the discrete entropy in 
a cell for any scheme, obtained in [22]. See equations (4.2), (4.3), and (4.4) 
below. In [22], a class of approximations called E schemes was shown to 
be convergent, even for nonconvex, but scalar flux functions. This seems to 
be the widest known class of convergent schemes in this general case, but 
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unfortunately the approximations are, at most, first order accurate. 

All the results in this paper are for semi -discrete (continuous in time) 
approximations, and thus can serve as guideline for a wide variety of time 
discretizations, both implicit and explicit. One of the principal applications 
of this theory has been in supersonic and transonic aerodynamics, where 
phenomena close to steady state are often studied. As an important example, 
one might mention the enormous advantage found when algorithms used to solve 
the transonic small disturbance equation, baaed on traditional Murman [2I] 

(Roe's unfixed for scalar equations) differencing were replaced by using the 
E schemes of Engquist and Osher [I2], [7]» later Godunov [13]. This 
minimal coding change, based on an entropy fix, increased the robustness 
of production codes by an order of magnitude, or more. 

Goodman and LeVeque, [11], have recently obtained a rather depressing 
result - two space dimensional, scalar approximations, cannot be TVE, and 
still be more than first order accurate, if the associated flux functions 
are reasonably smooth. Nevertheless, two dimensional schemes based on 
dimension by dimension high resolution differencing, do perform extremely 
well, even for complex configurations, with very strong shocks. See [4], 

[16], [3]> and Section VII below. This phenomenon remains to be justified 
analytically. 

The format of this paper is as follows. In Section I, we review the 
relevant theory of weak solutions for hyperbolic systems of conservation 
laws. In Section II, we do the same for the theory of approximate solutions, 
and rederive results about bounding the variation in Lemma (2 .I). We then obtain 
a maximum (minimum) principle and a nonlinear saturation result - Lemmas (2*2) 
and (2.3)* Tti Section III, we develop a systematic procedure used to obtain 
high resolution schemes, based on any three point first order accurate TVD 
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scheme. We use second order " upwlndlng," together with flux limiters to do 
this. In Section IV, we show how to obtain an entropy inequality, and hence 
convergence, which is our main result, Theorem (4.1). We use artificial 
compression or rarefaction operators (ACR). Harten introduced the notion 
of artificial compression in [I5]. Here we use the compression part to firm 
up shocks, and contacts, and the rarefaction part to enforce an entropy 
condition. Rather precise admissible bounds on the ACR term are given. 
Section V is concerned with existence of steady discrete shocks for certain 
high resolution schemes. In Section VI, we construct our high resolution 
scheme for systems, and then prove an entropy inequality: Theorem (6.2). 


Finally, Section VII gives numerical evidence demonstrating the utility 
of these schemes, as well as some results illustrating a few theoretical 
points made throughout this paper. See also [3]. 

Acknowledgement: The authors would like to thank Peter Sweby and Bram van Leer 

for some very valuable conversations concerning flux limiters. 

I. Review of Theory of Weak Solutions 

We shall consider numerical approximations to the initial value problem 
for nonlinear hyperbolic systems of conservation laws 

(l«l) fCw) =0, t > 0, -1 < X < 1, 

with periodic boundary conditions : 


w(x + l,t) = w(x,t). 
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and initial conditions w(x,o) = Wq(x). 

Here w(x,t) is an m-vector of unknowns, and the flux function, f(w), 
is vector valued, having m components. The system is hyperbolic when the 
Jacohlan matrix has real eigenvalues. 

It is well known that solutions of (l. l) may develop discontinuities 
in finite time, even when the initial data are smooth. Because of this, we 
seek a weak solution of (l.l), i.e., a hounded measureahle function w, such 
that for all 9 e Cq(R x R'*'), 


( 1 . 2 ) (a) 



f(w)9^)dxdt = 0 


(b) 


lim llw(x,t) - w (x)ll = 0 
t-^0 ^ L-^ 


Solutions of ( 1 . 2 ) are not necessarily unique. For physical reasons, 
the limit solution of the viscous equation, as viscosity tends to zero, 
is sought. In the scalar case, this solution must satisfy, for all 9 e Cq(R x R 
9 > 0, and all real constants c ; 

(1.3) (a) + sgn(w - c)(f(w) - f(c))9^)dxdt <0. 

This is equivalent to the statement: 

( 1 . 3 ) (13) ^ ^ ((f(w) - f(c))sgn(w - c)) < 0, 

in the sense of distributions. 

Such solutions are called entropy solutions. Kruzkov has shown*in [I 8 ], 
that two entropy solutions satisfy; - 

(l*^) lk(x,t^) - v(x,t^)ll ^ < l|w(x,tQ) - vCxjt^)!! 

Xi L 

for all t^ —^0’ condition (l. 3 ) guarahtees the uniqueness of solutions 
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to the scalar version of (1.2 ). Existence was also obtained in [18], 

For systems of equations. Lax has defined an entropy Inequality [19], 
with the help of an entropy function V(w) for (l.l), defined to have the 
following properties: 

(i) V satisfies 


(1.5) 


V f 


w w 



where F is some other function, called the entropy flux. 

(ii) V is a convex function of w. 

It follows from (l.l), upon multiplication by V^, that every smooth 
solution of (l.l) also satisfies; 

( 1 . 6 ) 

It was also shown in [19] > that if w is the bounded a. e. limit of 
solutions to the regularized equation, then the limit satisfies, in the weak 
sense, the following inequality; 

(1.7) \ + 

Inequality (l.3)(b), for scalar equations, is just (I. 7 ), with V(w) = 
w - c I . 

Ine'qualities (I. 3 ) (l. 7) have important geometric consequences for 

piecewise continuous solutions. Suppose w(x,t) is such a solution having a 

jump discontinuity, w^(t), w (t) moving with speed s(t). Then (1.2 ) 

h i\ 

implies the well known jump conditions 

(1.8) f(w^) - f(wj^) = s(w^ - w^). 

In the scalar’ case, (I. 3 ) is equivalent to Oleinik’s condition E across 


the shock 
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(1.9) 


f(w) - f(v^) 


w 


V, 


R 




Vt 


w. 


R 


for all w 'between w_ and w„. 

li K 

If f is convex, then (l.9) is equivalent to the statement that character- 
istics flow into the discontinuity, as t increases. Also, for scalar 
convex f, inequality (I. 7 ) for a single fixed convex Y is equivalent to 
( 1 . 3 ) for all constants c, if the solution is of bounded variation. This 
is a consequence of the recent results of DiPerna [ 5 ]. Thus uniqueness in 
this case is a consequence of Kf uzkov* s results [18]. This fact is crucial 
to our convergence proof in Section IV. 

For hyperbolic systems of equations, the Jacobian of f, denoted by df, 
has real eigenvalues, which are usually assumed to be distinct; 


^1 < ^2 * * * 


< X 


m 


> 


corresponding to right eigenvectors: r^,r^, . . . ,r^. Lax [19]> then defines 

the field to be linearly degenerate if 


V X, • 
w k 


= 0 . 


He also defines the field to be genuinely nonlinear if 


V X, .r, 4 0. 
w k k ' 

For genuinely nonlinear fields, a k shock moving with speed s, is 
defined to be a discontinuity of w, such that m + 1 characteristics 
flow into the shocks, and 

\("l> - > 0 > 

This geometric condition is equivalent to (l. 7 ) for weak shocks [19]. 
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II. I^ellmlnary Theory of Approximate Solutions 

We consider a semi -discrete, method of lines, approximation to (l.l). 
We hreak the interval [-1,1) into suhintervals 


'J. = {xl (J - ^)A < X < (j + |)a}, 

J 


J = 0,+l,...,+W, with (2W + l)A = 2. 

Ijet X. = JA, he the center of each interval J., 
J J 

Define the step function for each t > 0, as 


with end points 


V^(x,t) = Uj(t) 


for X € c9.. 

J 


The initial data is discretized via the averaging operator 


( 2 . 1 ) 


TaW-Cx) = i / w(s)ds = u.(o) for x e J.. 
A 0^ ' A ^ ' 0 <3 


For any step function, we define the difference operators : 


D u. = i A u.. 

+ 0 ^ + 0 

A method of lines, conservation form, discretization of (l.l), is a system 
of differential equations: 

(2.2) ^ u + D h 1 = 0, j = 0,+l,...,+N 

OO J + J ~2 — — 

^^(x,0) = T^Wq(x), for X e Jj. 

Here, the numerical flux defined by; 
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(2.3) 

for k > 1, is a Lipschitz continuous function of its arguments, satisfying 
the consistency condition: 


h(w,w,...,w) = f(w). 


It is well known that bounded a.e. limits, as A 0, of approximate 
solutions, converge to weak solutions of (l. l), i.e. (l.2)(a) is astisfied. 
However, this does not also imply that the limit solutions will satisfy 
any of the entropy conditions mentioned above, e.g. [8], [9]. Some restrictions 
on h are required. 

A simple class of flux functions h, for which (2.2) converges to the 
unique entropy solution in L°(L^(r); [0,T] ), as A -♦ 0, for any T > 0, 
are scalar "If schemes, [22]. Such schemes satisfy the following: 

Definition (2.I) . A consistent scheme whose numerical flux satisfies 

(2.4) sgn(u - u )[h 1 - f(u)] < 0, 

J J“2 ~ 

for all u between u. , and u., is said to be an E scheme. 

j-1 y 

It is easy to see that this class includes the widely known class of 
three point monotone schemes, i.e. those for which 


h T ) 

J J J 

with h. nonincreasing in its first argument, nondecreasing in its second. 
J"2 

Ib,rtial derivatives of a numerical flux, when needed, will be denoted 


via: 


8 


J+Y 
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Thus a three point scheme is monotone iff; 

\ ° < "^ 0 - 

One particular three point monotone scheme is due to Godunov [10], and 
has a special significance to this theory. The flux for Godunov* s scalar 
scheme can he defined hy 

(2.5) 1 ^ 

J “2 G 0 0-1 0-1 - 0 

0-1 0 

= max f(u), if u. ^ > u.. 

u. ^>u>u. ^ 

0 - 1 - - 0 

Thus, one can characterize E schemes as precisely those for which 

(2.6) (a) h. 1 < h^ if u _ < u 

(h ) h 1 > h^ I if u > u . 

It follows from [22], Lemma (2.I), that these approximations are at most 
first order accurate. 

Together with an entropy inequality, a key estimate involved in many 
convergence proofs, is a hound on the variation. We present an argument 
originally due to Sanders for monotone schemes, [29 ]> used to obtain this 
hound. 

For any fixed t > 0, the x variation of V^(x,t) is; 

j 


X. 1 = 1, if A u. > 0 
’ + 0 - 

X. 1 = -1, if A u„ < 0. 
0+i + J 


Let 
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Then: 


_d_ 

dt 


= X ) 

'' '' A dt ^ j+5 + j' 


1 y-x. , ;^AU. 

A Z-r j+f dt + j 


= _ J V X . lA A h . 1 
A J+2 + + J“2 


= J V (X. 1 - X. i)A h. 1 < 0, 
A ' j+t + J-i “ 


if we can write 


(2.7) (a) 
(c) 


A h. 1 = -[C. lA u. - D. lA a.] 
+ J-i J+2 + J 0-2-J 


C. 1 > 0 
J+2 - 


D. 1 > 0. 
J-2 - 


In the case of E (or 3 point monotone) schemes, this follows by- 
defining: 


C. 1 = 
J+2 


A u. 
^ + J 


D. 

J 


1 


f(u ) - h 1 

J J 


A u 


■J 


as in C22]. 

Harten in [I6], pointed out for explicit me-thods, that a variation bound 
could be obtained for schemes which' are higher order accurate. In our present 
method of lines context, it involves a five point consistent approximation 
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(2.8) 

(a) 


s=— — All. = lA u . — D . xA 1 

A + 0+2 + J " 

with 

(b) 



and 

(c) 


”j- 4 ' 


both Lipschitz continuous functions of their arguments. (See also van Leer 

[22]. ) 

We have thus proven for schemes of type (2.8). 


Lemma (2.I). For any t > t^ > 0 
' — — ^ ^ 

We also have the maximum and mi nimum principles : 

Lemma (2«2) . Let max u.(o) = M, min u.(o) = m. Then, for 0 < t 

J 0 

and each j : 

m < u (t) < M, 

t) 

Moreover, if C. 1 > 0, and u.(t) = M, then u. ,(t) = M. If D. i<0 
0+2 0 0+-L 0^ 

and u.(t) = m, then u. ,(t) = m. 

J 0 


Proof. The proof is trivial if the initial data is constant. Otherwise, 

iW 


we let satisfy 

(2.9) ^ = (C. 1 + e)A u^ - (D 1 + e)A u 


^ £ — fn _^A f T\ _^A 


for e > 0, with the initial data u?(o) = u.(o). 

0 0 

Suppose the maximum of u^(t), . for 0 < t < T, occurs for u? (t^^), 

0 ^0 ° 
with t > 0. Then (2.8), (2.9)? imply that u?(t^) = u? (t^). Hence 

0 \ /7 X X /J X- J ' 0 Oq 0 

u?(t) = constant, for 0 < t < T. This is a contradiction. Thus: 

t) 
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max 

j,0<t<P 




< M 


Let e 4, 0. It follows that u.(t) < M for each j, t. 

J 

The remainder of the proof follows analogously. 

Next, we prove a limit on the possible accuracy of approximations (2.8). 


Lemma (2«3) » Approximation (2.8) is at most first accurate at nonsonic 
critical points of u. 

(a sonic point u is one such that f’ (u) = 0. ) 
o 

Proof. For any C function u(x), we let C. i(u(x + 2A),u(x + A),u(x), 

J+2 

u(x - A)) be denoted by C. 1, C. i(u(x),u(x),u(x),u(x)) be denoted by 

J+2' 0+2 

C. i(u), and similarly for D. 1, 

0+2 J ~2 

Then 

(2.10) - ^A_h (u(x + 2A),u(x + A),u(x),u(x - A)) 

= (C., i.(u) - D._i(u))u + f u (C i(u) + D ^(u)) 

J+2 ^ ^ J+2 J“*2 

+ - C. i(u) - D 1 - D i(u)) 

X J+2 J+2 J~2 0~2 

+ o(a). 

Consistency implies: 

(2.11) C i(u) - D ,(u) = -f« (u), 

0+2 J "2 

while second order accuracy at critical points means 


(2.32) 


C. i(u) + D i(u) = 0. 
J+2 J~2 


Solving (2.11), (2.l2)j gives us: 
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J+2 

D, i(»i) = i f’ (li)- 

J"2 

Thus one of inequalities (2.8)(b), (c) must fail at nonsonic u. 

Thus, although schemes of the type (2.8) can he made to he as high as 
third order accurate, Lipschitz continuity of C 1, D. 1 implies a local 

tJ+2 0 "2 

degeneracy to first order accuracy at smooth maxima and minima. This local 
degeneracy, together with some results on initial boundary value problems in 
[l 4 ]. Indicate strongly that overall second order accuracy is the best 
possible. 


III. Total Variation Diminishing, Second Order Accurate, Scalar Approximations 


We now describe a systematic procedure used to construct second order 
methods of the form (2.8) from three point, first order methods of the same 
type. 


We begin with a three point scheme 


8u„ 


(3.1) ^ s 


= i c\^}a u. - i d4iA u. 
AJ+2+0 ^ J -2 - J 


- h(u ,u )) 
J J 


for 
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^ h(».^^,a.) - f(».) 

I A u. 

' + J 

f(.l) -h(.i.a^-^) v 

) 


> 0 


> 0 


These quantities are both nonnegative for E schemes. Moreover, the 

entropy inequality was shown to be valid in [22] for this class. Thus, using 

Lemmas (2.I), (2.2), and a standard argument, e.g. [291? 'w® have convergence 

to the unique entropy solution, as A 0; for this class. 

However, the nonnegativity of these functions, by itself, need not imply 

convergence, since the entropy conditions may fail, as it does for Murman’ s 

(Hoe’s) scheme. See e.g. [3^]. We do assume nonnegativity of c(^i, 

J+2 J "2 

in what follows, as well as Lipschitz continuity of h(u,v). 

Our first attempt at constructing a higher order accurate scheme comes 
from a simple upwind type of hybridization; 
h\i 

(3.2 ) 3? “ S - 5 * i 




with the numerical flux satisfying; 


(3.3) i = h(u ,u ) - ^ (h(u ,u ) - f(u )) + ^ (f(u ) - h(u ,u )) 
J “2 J J J J J 

(Throughout this work we shall often use the fact that consistency is equi- 
valent to the statement: 

h(u,u) = f(u). ) 


This scheme is, within second order accuracy, just the central difference 
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algorithm, and is thus, always second order accurate. It can not, by Lemma 

(2.3)j ot>ey (2.8). 

We can define the associated quantities so that 

-A u. - d(^)a u. 

+ J-2 J+2 + J J-2 - 0 


where, in this case: 


(3.4) 


0+5 V 


0+5 


1 ^^^“1+2’".1+1> - ‘^^".1+1’ "l+l^ . / 


D 


(2) _ Jl) 
J-4 ■ “o-i 



h(u. ^,u. t) -h(u. T,u. ^ 
^ j-r ,1-1^ ^ ,1-1’ J-2 

h(u., u. )~ - h(u.,u. T ) 

0 0 0 0-1 



It is cleax that the two quantities on the right above may become negative 
if the nontrivial ratios above become sufficiently large. 

One possible remedy is to employ a flux limiter, using the notation of 
Sweby [30]. (See also van Leer [32], for very original related work.) 

L<;t : 


(3.5) 


f(u ) - h(u ,u ) 

p _ J+-L ,1 J 

S - f('Oj.i) • 


Oirr TVD approximation will be of the form 
1 

(3*6) ^ [h(Uj^^,Uj) - h(Uj,Uj) - i 

This scheme is easily seen to be second order accurate away from critical 
points, (i.e. points where the denominators in (3.5) approach zero), if the 
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Lipschltz continuous function, satisfies ilt(l) = 1. 


The resulting quantities C. 2 ,, D. 1 are required to he positive: 


(3.7) c.i = -f J Li 


J+2 




A u. 
+ J 


h(.^,up - h(u,,u..,) 


A u. 

- J 


1 + 


1 + i 


rt(RT) 

“ - I I > 0 


L R. 

- Rt 


«>]] 


- ^ 0 * 


We thus have the required inequalities for '1 i(R— ): 


( 3 . 8 ) (a) 

M 


1 + ^ 


1 + i 


KrT) 

u 

L r7 

J 

KRp 

J 

L Rt 




^(R'l) 


> 0 


> 0. 


Various slope limiters have been developed. See Sweby [ 30 ], for a 
numerical and theoretical analysis of their properties. As an example, we 
may take 


( 3 . 9 ) 


KR) =0 if R < 0 

ilr(R) s R, if 0 < R < k for 1 < k < 2 
^(r) = k, if R > k. 


Let the numerical flux defined from (3«l)“(3*8) be called 




Its precise definition is 


(3.10) - i - h(Uj_^,Uj_3^)) 
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We thu£5 have the following: 


Theorem (3»l) . 
( 3 . 11 ) 


The approximation defined through 

^u. , 




WQ(s)ds, 


with Wq € n L° n BV, defines u^(x,t) having the property, that as 
A 0, there exists a converging suhseq.uenee (x,t) which converges in 

L (L'^(R), [0,,T] ), to a weak solution of (l. l). 


The proof is a routine consequence of Lemmas (2.I) and (2.2). See, 

e.g. [29]. 


Theorem (3«2) . The approximation (3.I0) is second order accurate, 
except at isolated zeros of h^ (w,w)w , or h (w,w)w . 

X D iK 
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IV. High Resolution Schemes and the Entropy Condition for Scalar Approximations 


Following recent tradition, we christen the schemes constructed in the 
last section high resolution schemes - namely they are second order accurate 
(with the usual exceptional points), variation diminishing, and have a five 
point bandwidth. The question remains - is it true that limit solutions 
satisfy the entropy condition? For first order non E schemes, variation 
diminishing is not enough. The perennial example is Murman' s (Eoe’ s ) scheme 
with f"(u) >0 and 

(4.1) h^(Uj^^,Uj) = I (f(u^^^) + f(u-j)) - i 

L R L R 

Given an entropy violating shock u < u , with f(u ) = f(u ), it is 
well known that 


A f(u.) 

+ il 


A u. 
+ J 


A u. 
+ J 


u. = u 


for j < 0, u. 


R 


j > 0 


L R 

is a steady solution to (3.I). Moreover, since h(u. ,,u.) = f(u ) = f(u ) = h(u.,u.), 

J+-L J J J 

in this case, then any of our high resolution schemes (3.II), will have entropy 

M 

violating solutions, if h(u. ,,u.) = h (u. ,,u.). 

J J+J- J 

Let V(w) be any convex function. In [22], Section III, it was shewn 
for any solution of (2.2), that 


(4.2) A(— V(u^) + D^F^(u.)) 


y '* J+1 
u. 


u. 

J 


dwV" (w)[h 1 

J+2 


- f(w)l. 


where the approximate entropy flux is defined through: 
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ih.3) 


F (u ) = F(u ) + V'(u )[h - f(u )]. 

J J J 2 J 


Thus, a sufficient condition that any limit solution satisfy (1.6), for 
a fixed V, is that 




■j+l 


f V”(w)[h 1 - f(w)]dw < 0. 

J \ X . '^•*■2 


In order that the above Inequality be valid for all convex V, it is 

necesssiry and sufficient that h. i correspond to an E scheme, and hence 

J+2 

that the approximation be at most first order accurate, [22]. Thus we shall 
only obtain our entropy inequality for a single V(w), say V(w) = -g- w^. 

This is sufficient for convergence, if f is convex, as will be shown in 
Section IV. 

We now proceed to modify the flux difference quantities, to ensure that 
is valid, and, moreover, to sharpen discrete shock profiles. We shall 
do this by usin^ the notion of artificial compression introduced by Harten 
[ 15 ], cf also [ 16 ]. We also add negative artificial compression (= artificial 
rarefaction), if f’ (u.) < f* (u. ^ ), and a certain amount of (positive) 
artificial compression if f’ (u.) > f’ (u. ^ ). The hi^ resolution properties 
are preserved, and (4.4) is shown to be valid. 

It is known that too much artificial compression can cause expansion 
shocks to develop. We shall obtain fairly precise bounds on the amount allowed. 
Once again, we define 


^(iij'v) = h(u,v) 


h(u,v). 


We shall also use two Lipschitz continuous functions of u. u., 

d+1 d 
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*1* 

denoted by a. i, a. i and defined below. Our scheme is an artificially 
•3+2 3~2 

compressed or rarefied version of (3*6). 

We first modify our flux differences via the following: 


m 


(4.5) (a) (h(Uj^^,Uj) - h(u^,u^)> 

tlrJ- J J J 


1 + a. 1 

- “ njLl. 

' ^ J+1’ J 


J+2 




h (\i ,VL ))A \1 

i_x_4 

h(u,,u,)) 


J J 


ih) 






.4. 

Here a. 1 are both positive numbers, chosen first so the quantities in 
3+2 


brackets in 

(k.5) (a,b) are both always between 

0 and 2, i.e. we must 

have; 





(4.6) 

(a) 


VI 

1 

(h(u ,u ) - h(u ,u )) 
J+-*- 0 J J 

A u. 

+ J 


(^) 





A u. 
+ J 


We now let; 


(h(u,,up -h(u.,u^-l)r 

' (‘“'Vi’Vi' ■ 

m (h(u ,u ) - h(u ,u ))“ 

R " = Jt-I- . 3 J J 

^ (h(u.,u. t) - h(u. T,u. ,))™ 

' ' j’ 0-1^ j-1’ 0-1 


Our scheme, which uses 


ACR (artificial compression-rarefaction) is: 
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Sa. m 

(4.7) [h(u.^^,u.) -h(a.,u.) -iMKY-L)(h(u.^l,u.) -h(u.,u.)) ) 

m 

+ ) - h(u ,u A^((i1i(R ■J^^)(h(u ,u ) - h(Uj,u .j^)) )] 


where we take ')i(E) to he as defined in (3.4), wi.th k = 1. 

We let the numerical flux defined through (4. 5)j (4.6), (4.7) he denoted 


hy H^°i. It is precisely defined by (3.IO), with the superscript m 
<J“2' 

attached to both Rv and the two flux differences. 


Now, for TVD, we must have: 
h(u, ,,uj - h(u^,u, ) 


0 < - 


j+1’ .1 


A u. 
+ J 




(h(u^^l,Uj) - h(u,,u, )) 


m 


m 


1 i J+-J-' J • ,r .1 

^ ^^( - h(u.,uj) 
J+-L J J J 


h(u ,u ) - h(u ,u ) 

0 < ( J . V j . 

A U. 

- 0 


(h(u ,u ) - h(u.,u._T ) 

V 

.. T ,U. 

J+1 J 


m 


^ Ih(u'^ ,^uj - h^u,,u, )) 


r1;(R ") m 

m j+1' 

R." 


in 

4-^ - ^(R.%) 

m 0-1' 

R + 


which is valid, because h satisfies (2.7), and because of (4.5), (4.6) 

and (3.10). Ttieorems (3.I) and (3.2) are also valid for this ACR version of 

our scheme, but our real goal is to choose the av 1 so that the inequality 

0 + 2 


(4.8) 


r 0+1 
J \i. 


- f(w)]dw < 0 
0 + 2 ' ~ 


is valid, i.e. to enforce inequality (4.4), with V(w) = k 


w 


It is easy to see that 


ih.9) 



-i 

*^0+1 

+ t(^(A^Uj))^ - (w - + up)^]f"(w)dw. 

^0 
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Substituting this into (4.8), gives us the equivalent, desired inequality 

(4.10) i J [(i A^u^^ - (w - + Uj))^]f'(w)dw 

m 

- b4i'h.(Vi’V2> - 

m 

+ i(A^u.)[Ch(u^^^,Uj) - h(u.^^,u.^^))“(l - ♦(R^*)) 

From (4.5), (4.6), and the fact that h is an E flux, it suffices to 
prove (after cancelling the factor h) 


Ml) / 

n 


f" (w)dw[(i A u.)^ 

+ J 


(w - 






so- 


We thus have; 


Lemma (4.l) . Any solution to (l.l) which is the limit of a subsequence 
of approximate solutions Ir^, as A' -♦ 0 of 


(4.32) 




WQ(s)ds, 


satisfies the entropy inequality (l. 6), provided that inequalities (4.6) 
and ( 4 . 11 ) are valid. 
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As mentioned in the introduction, for scalar convex conservation laws 
where solutions lie in BV, a single entropy inequality implies uniqueness 
of solutions. This follows from the recent results of DiPerna [5]. We know 
that if w(x,0) € BV, then the function h^(x,t) has x variation hounded 
hy that of w(x, O), for any fixed t. In particular 


+ h,t) - V^(x,t)ldx < var(w(x, 0) ) 


for any t > 0, A > 0, h > 0. 

Tlie dominated convergence theorem then guarantees that any limit solution 
w(x,t) will have x variation hounded hy w(x, O), uniformly in t. Consider 
the interval 0 < t < T, for any T > 0. Then 
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Thus the dominated convergence theorem implies that any limit solution 
w(x,t) will have bounded variation in x and t. This together with Lemma 
(i|-.l) and the above mentioned uniqueness result gives our 

Theorem (4.1) . (Convergence) The sequence of approximate solutions 
converges a.e. as A -► 0 to the unique solution of the scalar convex 
conservation law (l. l) provided that the initial data is in BV and that 
inequalities (4.6) and (4.11) are valid. 


For simplicity, we exemplify this theory by considering semi -discrete 
versions of the following three monotone schemes. 

(4.15) (a) Lax -Friedrichs ; 

‘'“''Vrt’ “ * i 

with 0 < k chosen so that |f* (u)| < k6, for some 0 with 

0 < 0 < 1 . 

(b) Engquist -Osher : 




Eo /* j+1 r j 

h (u. T,u.) = / min(f* (s),0)ds + / max(f’ (s),0)ds - f(o) 

^ Jq Jo 




(c) Godunov: 


G 


h (u. ,,u.) = min f(u) if u. < u 


j+1’ 0 


u.<u<ii. T 
J J+1 


■j j+1 


max f(u) if u. > u. 


u.>u>u. T 
J 0+1 


■j 0+1 


For lax -Friedrichs, (4.6)(a) and (b) become 
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(4. l6) 


'Vl> ■ i 

J 


while (4.1l) is valid if; 


(^.17) 


r/ A ] 
(w )dw ( —4 


1 ) -(w-*(u.^,.u^))^; 


< min n't i(A u )2(ft (u ) - f (u )). 
“ . 3+z + 0 <J+-i- J 

t-« “ 


Cleairly these inequalities are compatible if 6 is sufficiently small, 
fforeover, if f is convex, then (4.16) becomes 


(4.18) (a) if Uj < (rarefaction) 


(w)dw ( 


(w - -Ku. + u. ))^ 
J+1 0 


at , > 
j+i - 




(b) if u. > u. T (shock) 

J J+1 


at 1 < — i 
0+2 - 


r'^O+l r / 

J f"(w)dw[(^-t 


(» - H. 


<Vj) 


The right side of (4.18) is always positive for convex f. For example, 
for f(w) = -I w^, the right side is ^ . 

The precise restrictions for this case become 


(4.19) Burgers’ equation 


laj+iVjl “j' 


and 
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Thus the high resolution scheme based on LF is a convergent approximation 

i. 

to Burgers' equation if 6 is small enough and the then chosen as 

above. 


For the Engquist-Osher scheme, (4,6)(a) and (b) become; 


(4.20) 




(f (u, ,) - f (u,)) 


j+2 - J+1 


while (4.11) is valid if: 


^.21) f 

^ n 


j+i r / 

f»*(w)dw -±-i 
+ _ \ 2 


(» - + aj))2 


J+2 j + <]+■*- “ J 


Inequalities (4.20) and (4.2l) are obviously compatible for convex f, 
for certain bounded a~ i, if !a u.l is sufficiently small (i.e. in a 

<3+2 + J 

region of smoothness). They are also compatible for 1 a u.j bounded away 

+ J 

from zero. This follows since we may take in (4.20) 


(4.22) 


+ 

aT 1 = 
J+2 


qr (A f u.) 

^ ^ + T- ,1 [ 

(A_^u. )(aV (u 

-r J T -f. J 


It now remains to prove that (4.2l) is valid, i.e. 


/ Vl r/Au.v^ -1 

f"(w)dw^-!^ j - (w - + up)2 + (A^f^(u^))A^u^ < 0 


or equivalently 
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(4.24) 


/ J+1 

- f_(w))dw < 


u. T 

r 0+1 

/ (f (u.) - f (w))dw < 0. 

A. + J 


These last are obvious facts. 

In the case f(w) = -g then (4.2l) becomes 

(4.25) (Burgers’ equation) 

(a) if Uj < (rarefaction), then: 


(i) 


a? 1 > I 

0+i “ 6 


at 1 > 0 


0+2 


if u. ^ < 0 
0+1 - 


(ii) 


a 1 > 0 
J+2 

a . \ ^ v 
0+4 “ b 


if u. > 0 
0 - 


(iii) 


a” , > i _ 1. f ^^0+1 ^0 \ 

M-4 :12 I )2 ; 

+ 0 

? 2 

1 1 / 3u.. - + u. 

at , > 1 + ,1- ( _J±1 si 

0+2 12 32 \ j2 

+ 0 


if u. , > 0 > u. 
0+1 0 


(b) if Uj > (shock), then 


( 1 ) 


^0+4 - 6 


at 1 < 0 

0+4 - 


if u. < 0 
0 - 


(ii) 


a. 1 < 0 

0+4 - 


+ 

^0+4 - 6 


if u. T > 0 
0+1 - 
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- .1 1 


12 


(iii) 


3u^ + 

J J+1 

(A u.)^ 

' + J 


+ .11 

j+i - “ 12 32 


2 2 

/ 3u. + u. T \ 

V (au.)2 ; 

+ 3 


if u. T < 0 < u.. 
J+1 3 


Moreover, if jat ij <-5- (which is admissible from (4.25)? then (4.20) 
3+2 - 

is no restriction at all. 


Thus, the E-0 scheme incorporated into H. i yields a convergent high 

3+2 

resolution scheme for any convex f, and any initial datum having bounded 

+ 

variation, provided the aT are judiciously chosen. 

J+2 

A tedious calculation yields the analogous results for Godunov’s scheme. 


V. Steady Discrete Shocks 

We now check for the existence of discrete, steady, shock solutions to 
the high resolution, entropy condition satisfying schemes constructed in the 
previous section, based on the Engq.uist-Osher flux. Since the scheme satisfies 
Lemma (2.I), the profiles must be monotone. We shall also show that they are 
sharp. 

It R 

Let u , u , be the left and right states for a physically correct 
steady shock, i.e. 

f(u^) = f(u^), 

Ii I» R 

f(u) < f(u ) for all u between u and u , 


A steady, discrete shock {u.} will satisfy: 


lim u. = u 

j-4U» J 


11m u. = u 

j_-oo J 


R 


(5.1) 



31 


for all j. 

0+2 


For simplicity, we assume f"(u) >0, with f’ (o) = 0 and f(0) 
The numerical flux becomes 


= 0 . 


M 


(5.2) ■ 4 


M 




If we take k = 1 in (3.9)? we have 


M 

( 5 . 3 ) (a) 


M 


-•= u max[0, min ( [A f (u. + a“ i _ ^(A f (u )A u ]( ^ u )'^) ] 


+ 0 


r=0,l 


M 




H - 0+^ O+t+3^ + - 0+^ + 0+3^ + 0 


M 


A u. max[0, min ( [A f (u. 
+ r=0,l 




We shall seek steady discrete shocks of the same general form obtained in 
[8] for the E-0 scheme. 

These are 

(5.3) u = u^ j < -1 

t) 

= U^, j > 2 

> Uq > 0 > > u^. 

For the first order scheme, u. can be viewed as a smooth function of 


u^, satisfying the above and 


f(aQ) + f(u^) = f(u ). 
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For the present scheme, we shall get a different one parameter family of 
intermediate states. 

It follows for all j ^ 0, that (5*3) implies; 

For j = 0, we have the following equation; 


(5.4) 0 = f(u^) + f(uQ) - f(u^) 

f(u^) +alf» Uq) f(u^)-f(u^)+a" (f (u^)-f' (u^))(u^-u^) 




R, 




\ 


0 


-f(up) - ati(-f»(uQ))(u^ - Uq) 
^ - ^0 


= G(u^,Up). 


R 


One special solution is, again, u^ = u , Up = Q. We have 

(5.5) 




% (»^ 0 ) . 0 




(u^,0) = ( ^ - i a" )f"(u^) > 0 
au| 2 3/2 

, R „x 


^G _ i 


2 = i f"(0) > 0. 


A simple application of the implicit function theorem gives us; 


Theorem ( 5 . 1 ) . There exists a family of sharp, discrete, shock solutions 
to ( 4 . 12 ) of the type ( 5 . 3 )? with u^ a smooth function of Up and 0 < Up 
small enough. By symmetry, the same is true with Up a smooth function of 
u^, and -u^ small enough. 
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VI. Systems of Conserva.tioQ Lavs 


We shall now build a second order accurate scheme approximating (l.l) 
for systems: m > 1 and verify some of the desirable properties mentioned 
above. Our basic three point first order scheme will be the one devised 
by the first author [2^], with Solomon [26], and analyzed jointly by us in 

[ 2 ], [25]. 

The numerical flux for this scheme is constructed as follows. First, 
we define a piecewise smooth, continuous path in phase space, connecting 
w^ to made up of m subpaths. Along subpath k, we have 


( 6 . 1 ) 




*ttl 

with r, , the k eigenvalue of 6f(w) corresponding to eigenvalue \ . 


■k’ 


(Here \. < X < • • • < X . ) 
' 1 2 m ^ 
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On such a subpath, the m - 1 independent Riemann invariants, corresponding 
to field k, are constants. Our ordering is as follows. ¥e begin at w 
with k = m, stop at some endpoint, use k = m - 1, etc., arriving at 
w. ^ with the subpath corresponding to k = 1. ¥e call each such subpath 

J + l n 

V V 

rl 1, Y = m,m - 1 , ..., 1 . There exists exactly one such path R. i = U Ri i.> 

j+2 Y=m 

for j 1 sufficiently small (actually, for physical systems such 

as Euler’ s equations, the construction works in the large, as long as 

cavitation does not take place [26]). 

¥e define the numerical flux to be 

(6.2) h° 1 = I - f l^(w)jdw . 

J+2 J+-^ . J fr 

j+i 

Here ] ] denotes the absolute value of a diagonalizable matrix. 

This algorithm yields a closed form expression for h? 1 for various 

J+2 

physical systems, including Euler' s equations for compressible, Inviscld, 
gas dynamics, [26], [3], because their Riemann invariants can be easily 
tabulated, and because all fields are either linearly degenerate, or genuinely 
nonlinear. ¥e assume this property throughout the remainder of this paper, 
for (1.1), 

On such a k subpath, the integral in (6.2) is given by a simple 

expression involving +f(w) at endpoints of the subinterval, and, in the 

genuinely nonlinear case, it may also involve +2f(w), where w is the 

unique sonic point for which X, (w) = 0. 

K 
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The full algorithm can be interpreted, geometrically, as solving the 
incoming Riemann problem in phase space, using only rarefaction, compression, 
or contact waves, then averaging the resulting multivalued solution, in 
physical space, as in Godunov’s method [34], 

The first order scheme may be written so as to yield a field by field 
decomposition; 


( 6 . 3 ) 


■with 




(6.4) (a) h°(w ,w ) - h^(w ,w ) = 

0 0 0 


0 } 


X 


(6f (w) )'^dw 




(b) h°(w w ) - h°(w ,w ) 
J J J J 


y (6f(w))“dw. 
J 2 


Here 

(Sf)* = i ((Sf) + larl) 
(ar)- = i ((at) - iSfi). 


Fo 37 a linearly degenerate k field, we have X constant on R, xj 

k 3+2 

and hence: 


( 6 . 5 ) 


f (af)+au » 

✓ _k 




f (^)"dw = (1 - X(\^+2))[f(w'") - f(w^)] 

./■pk 


where X(x) =1 if x > 0, X(x) = 0 if x < 0, and w*^, w^ are the upper 


and lower limits of Integration. 
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( 6 . 8 ) 


5w 


1 

A 








k,j+i 


k=m 


fc=m 


+ ^ A 




r+ 


k, j' 


This differencing, although formally second order accurate, again leads 
to overshoot, undershoot, and occasionally nonlinear instability. 

To remedy this, we limit the flux differencing, as in the scalar case, 
except that we do this separately in each field. We shall use the entropy 
function, discussed in the first section, in the construction of the scheme. 
We remark, that although the existence of a convex V(w) follows only if 
the overdetermined system 


has a convex solution, the fact is that most of the equations of physics are 
endowed wiLth such an entropy. See [6] for further details. 

The entropy gradient will be used to construct a linear functional 
applied to vectors w. This is correct dimensionally, as pointed out by 
Harten-Lax in [I7]. 

We first define 

(6.9) V (w. , , 0 - V (w. , , \) = A^j’^V . 

' w' j+l-^k-l/m)"^ w' j+l-(k/m)'^ + w 


The algorithm in each field will involve the quantities 


( 6 . 10 ) 






.T+ .(j^k). 




X ) 


w' 


_ k, j+p + w 


> J+p 
Is 


(W+ i.a"'^^ 

K, J-2 


V ) 


for each field k = 1, ...,m. 
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Here ( , ) denotes the usual inner product in R™. 

Using any flux limiter of the type described in Section III (i.e. iJ;(r) 
satisfies ( 3 * 8 )> with Kl) = d)» w'e define a second order accurate, TVD 
type discretization: 


= - A ( E - i 

k=l 


m 




k=l 


We now show that for linear systems, this nonlinear scheme decouples 
to a second order accurate variation diminishing algorithm, and hence 
approximate solutions converge, as A -♦ 0, to the unique solution of (l. l). 
(See also [l6], for the first proof of such a result. ) 

We have 

df(w) = A = RAIi, 

with: 


and 


R = [r^ . . . ,r ] 

V 2 m-' 



the constant matrices made up of left and right eigenvalues of A, suitably 
normalized so that 

) = 0 if 1 ^ j 

i J 

= 1 if i = j. 
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and 



Then 


j )inax(X^, 0 

Moreover, since the entropy is just defined by V(w) = ■§• (wjw) for 
linear systems, then 


+ w 


(<ej^,A^w)r 


k* 


Thus, if we take the inner product of (6.11) with and denote 

(k) 


ty 


w; 


we obtain 


^wW = - I [min(X^,0)(A^w^’^) " i 


,(k) 


(k). 


+ max(Xj^,0)(A_w^ ^ 


This scheme is of the type (3*6)-(3«8)> and hence is TVD. Thus we have 

(k) 

the variation bound for each w. as in Section II, as well as the maximum 

0 

(minimum) principles stated there. ¥e thus have the simple; 


Theorem (6.1) . The approximation to the linear hyperbolic system (l.l), 
defined by (6.11), with initial data in n L° n BV, converges as A -* 0 
in L°°(L^(.'R), [0,T] ), to the unique solution of (l.l). 

Next, in order to sharpen shocks and contacts, and to enforce a discrete 
entropy condition for systems, ^we use ACR, as in the scalar case. We again 
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(as in the scalar case) require that iIi(r) be defined by (3-^)» with k = 1. 

In order to add ACE, we use the k eigenvalue difference along a 
k subpath. This is defined by 

(6.32) 


Next we define 


M 


(6.13) (a) N+ 1 

■K, J+2 


= n; 


k, j+-| 




^ + ’ + v' 


kfj+s'' + /JJ+ A(j7k)y ) 

' k,j+^’ + w' 


for genuinely nonlinear k fields, and 
M 




for linearly degenerate fields. 
Analogously, we define 


^(j,k).^ ) 


M / (A'— w,A--- 

(6.14) (a) a a;j^.(Ap> >X.) -t ^ 

' Ir n I i_ 


k,j+-5’ + ^w^ 


for genuinely nonlinear fields, and 

(t) 1 = N" id + a" . )^) 

k, j+2 k,j+2 k, j+2 + + w 


for linearly degenerate fields. 

We make the following restriction on the Lipschitz continuous functions 

%0+i' 

Restriction (6.1). The values of the coefficients multiplying wf . i 

k, j+2 

in (6.13) and (6. l4) are always strictly between 0 and 2. 
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M + 

We now define as in ( 6 . 10 ), replacing IC . i_ these compressed 

K, J t;, j+2 

or rarified flux differences. The TVD and ACR scheme is constructed as follows 


m 


„ MM 


k^l 


m 


k?=l 


1 A 

- ^ A 3+ . I • 
A - j+i 


We shall now prove our last result. 


Theorem (6.2). There exists values of ar . i , for each . 1 , k, so 

k, J+g- 

that if the solutions to (6. I5) with the usual initial data have, for A 
small enough, (a) sufficiently small oscillation, and (h) the property that 
eigenvalues corresponding to linearly degenerate fields stay hounded 

away from zero, then hounded a.e. limits of the approximate solution satisfy 
the entropy Inequality (I.7). 


li’oof . It was shown in [22] that it suffices to obtain the inequality 
(6.16) f (< 3 .w)\^(w)[M^°^ - f(w)] < 0, 

for a].l J, under the given hypotheses. 

An integration hy parts and rearranganent of terms gives us an equlvaJ.ent 


problem; 
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( 6 . 17 ) f (dw)^(^(w))[V^(w) - i + V^(wp)] 


- i 5: 




- i E 


+ i f(A V (w.), V' (k“ . 1 - N ". i)| < 0. 
^ V + ^ J Z_/ ^ k,j+-| k,j+^7 - 

' 1_ -1 ' 


Let k’ denote indices corresponding to genuinely nonlinear fields, and 
k" denote the remaining fields (which are all linearly degenerate). 


Then 




(6.18) (a) V(^ ^ ,i.\ 

/•u\ 1,± _ „i /.(j,k’7 A(j,k'7„ 

+ 

(Notice that the aT„ . x are not dimensionless quantities - they 

k , j+2 

— *4" 

have dimension inverse to that of (v)^. The aT, . i are normalized to 
he dimensionless. ) 

Using the argument (and notation) in [22], Section III, it is easy to 
show that 

(6.19) i (Wj), (^(w))dwj 

= i \ J <3.t(rj^(w(t))'^V^l\j^(w(t))lrj^(w(t)) 

^ 0 

m s 

+ i X) / / ^•+'(V^('^(s))2’^('w(s)) - V^^(w(t)r^(w(t)))^ 

v=i 0 Jo 

• l>^(w(t))lr^(w(t))dt. 
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let e. = w.l be sufficiently small. Then our assumption about the 

j + t) 

fields, the argument in [22], together with (6.20) and Restriction (6.1) 
guarantees that the second and third terms in (6. I7) are non -positive. More- 
over, regarding the last two terms, we have, using (6.I8), and letting 
C > 0 be a universal constant: 


/ “ M “ M \ 

(6.20) (a^V"3» £ <’‘^,0.4 - \:j-4> " S - \j.47 

' fc=l ksl 


k* 

. >V )] 

+ ’ + v'J 


k’ 




1)) 

2 


(A^-^’^’ V )2)|, 

' + ’ 4 - w / 


recalling again, that the k’ are the genuinely nonlinear, and the k" are 
the linearly degenerate, fields. 

Here we have taken 

(6.21) (a) a^, , i(Ap’^’ \ ) > 0 

y J'^2 i. 

(b ) aS, 1 < 0. 

» J+2 

This mesins that we have added only expansion, not compression, to our 
scheme. The freedom to add compression at shocks or contacts is lacking, 
only because of technical points in our proof. 

The Integral in (6.I7) can be rewritten as: 
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( 6 . 22 ) 


((dv)-*-V^(v))^(w)dw) 


. 1 ' • 1 TT 

0+2 0+2,^,U 


if if 

•'r,, 1 V •'r,. 


((dv) V^(v))^(w)dw)j. 


Here, F . i 


<3+2 


iicxc, X , 1 tt is that part of P. i which starts at w and ends at 

0+i,w,U j+f 

w. r. 1 T starts at w, and ends at w. Thus T. i = F i , - tt* 

J+1 0+i>'^}^ 0+2 t+2>^i^ 

Using the arguments of [22], Section III, we can show that the sura of all 


the integrals along subsegments is 




f Uj^(w(s))lds ^ , 


except, perhaps for a term: 


Sj^ / s 

(6.23) i i: /o <K - -f <it^( (rj^(w(t ) ) )'^V^(w(t ) )r^(w(t ) )\^^(w(t ) ) ). 


Let t = Sj^ - t* in the second integral. Using our hypotheses about 
the fields, easily gives us an upper bound 


m Sj^ 


for this integral, where the D are certain universal positive constants. 

k 

Comparing (6.23), (6.2^) to the right side of (6.20) shows that we may 

choose each constant a~ . i subject to Restriction (6.1), and sufficiently 

K, j+2 

large in magnitude, (if e. is sufficiently small) so that ( 6 .I 7 ) is valid. 

0 
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VII. Results of Humerical Experiments 

Since our preceding theory concerned only semidiscrete, continuous in 
time, approximations, we "begin by describing the time discretizations used 
for all calculations presented here. That was a two step time differencing 
found in Richtmyer and Morton [27]* In one space dimension it can be written 
as 

(7.1) (a) w’^'*'2 = ^ {first order approximation to f*^} 

2 X 

(b) = w’^ - At {first order approximation to 

- — {"upwind" approximation to +f*^ using flux limiters}. 

The resulting explicit algorithm is not, in general, truly variation 
diminishing, but still works well, even in multidimensional calculations. 

We begin by computing a solution to inviscid Burgers' equation which is 
a shock moving with speed The results are shown in Figure 1 . Next we 

compute a solution to the steady state Burgers' equation having an inhomogeneous 
right hand side, and compare the result with the known exact solution - Figure 2. 

Both of these scalar experiments were done using the flux limiter in 
( 3 » 9 )j with k := 1 , and with the high resolution scheme based on Engquist- 
Osher first order differencing. 

Next we compute solutions to Euler' s equations of compressible gas 
dynamics, using the schemes constructed in the previous section. See [ 2 ], 

[3], for the precise equations, and algorithm used here. Figure 3 shows 
the results of quasi -one -dimensional Laval nozzle flow. Reference [2] has 
the equation with precise source term. 
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We next consider the case of an oblique shock reflecting from a flat 
plate, in two space dimensions. Figure 4a shows a cartesian grid with 
6l X 2l points. Figure 4b shows a first order accurate solution, 4c was 
obtained using a non-TVD second order extension, and 4d shows the contours 
obtained using the TVD second order accurate scheme. Figure 4b has highly 
smeared incident, and reflected shocks, 4c has thick profiles because of 
undershoots and overshoots near the shock, and 4d has relatively sharp, 
nonoscillatory profiles. 

In Figure we present the computational grid to be used for a blast 
wave problem solved as an unsteady problem. In Figure 6, we give the computed 
pressure contours using the TVD, second order scheme, at time T = 0.469. 

In Figure 7, we show a finer computational mesh, and in Figure 8, we present 
the pressure contours using the TVD, second order scheme at time T = 0.459* 
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